home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 12 / Mac Magazin and MacEasy Magazine CD - Issue 12.iso / Sharewarebibliothek / Anwendungen / Wissenschaft & Technik / Yorick / yorick11-fpu folder / include / demo2.i < prev    next >
Text File  |  1995-06-29  |  8KB  |  270 lines
  1. /*
  2.    DEMO2.I
  3.    Mesh plotting demo.
  4.    To run, start yorick, then type the following two lines:
  5.  
  6.      #include "demo2.i"
  7.      demo2
  8.  
  9.    You can rerun the movies as many times as you wan by typing "demo2".
  10.    To run only the first, second, or third movie, type "demo2, n"; for
  11.    example to run only the third movie, type
  12.  
  13.      demo2, 3
  14.  
  15.    $Id$
  16.  */
  17. /*    Copyright (c) 1994.  The Regents of the University of California.
  18.                     All rights reserved.  */
  19.  
  20. func demo2(which)
  21. /* DOCUMENT demo2
  22.      Exhibit quadrilateral mesh plots in 3 movies of a drumhead.
  23.      The drumhead is initially stationary, but has a bump near one
  24.      edge.  Yorick is solving a 2D wave equation to compute the
  25.      evolution of this bump.
  26.  
  27.      The first movie is a filled mesh plot with color "proportional"
  28.      to the height of the surface of the drum.  A few well chosen
  29.      contour levels (here 3) add a lot to a filled mesh plot.
  30.  
  31.      The second movie is a "3D" perspective plot of the height of the
  32.      drumhead.  In this movie, the mesh lines are drawn, which is
  33.      slightly confusing since the cells are not all the same shape.
  34.  
  35.      The second movie is a "3D" shaded plot of the height of the
  36.      drumhead.  Yorick computes surface shading based on the angle
  37.      of each cell from a light source.
  38.  
  39.      As you watch this, you might reflect on the two dimensionality
  40.      of your retina.  What Yorick lacks by way of 3D graphics is
  41.      really just fancy hidden surface algorithms; the simple
  42.      painter's algorithm used here and in plwf.i is easy to
  43.      implement.
  44.  */
  45. {
  46.   npass= 200;
  47.   /* generate a 30-by-30 cell mesh on the [-1,1] square */
  48.   x= span(-1, 1, 31)(,-:1:31);
  49.   y= transpose(x);
  50.   /* map the square mesh into a mesh on the unit circle
  51.      this mesh has more nearly equal area cells than a polar
  52.      coordinate circle */
  53.   scale= max(abs(y),abs(x))/(abs(y,x)+1.e-30);
  54.   /* note that abs(y,x)=sqrt(x^2+y^2) */
  55.   x*= scale;
  56.   y*= scale;
  57.  
  58.   f= f0= exp(-8.*abs(y+.67,x+.25)^2)*(1.-abs(y,x)^2);
  59.   fdot= 0.0*f(2:-1,2:-1);
  60.  
  61.   af= abs(f(2:-1,2:-1));
  62.   lf= laplacian(f, y,x);
  63.   dt= 0.08*sqrt(2.*max(af)/max(abs(lf)));
  64.  
  65.   window, 0, wait=1, style="nobox.gs";
  66.   palette, "heat.gp";
  67.   limits;
  68.  
  69.   fma;
  70.   plt, "Yorick         (c) 1994 Regents of the U.C.\n"+
  71.     "X Division, Lawrence Livermore Lab\n"+
  72.       "\nInterpreter for scientific computing\n"+
  73.     "UNIX and MacIntosh versions\n"+
  74.       "anonymous FTP from:\n"+
  75.         "ftp-icf.llnl.gov:/pub/Yorick\n"+
  76.           "\nThis demo shows three movies\n"+
  77.         "of a drumhead oscillating",
  78.         0.40, 0.64, justify="CH", font="helvetica", height=24;
  79.   redraw;
  80.   pause, 8000;
  81.   fma;
  82.  
  83.   /* roll the filled mesh movie */
  84.   if (which && which!=1) goto persp;
  85.   fc= f(zcen,zcen);
  86.   cmin= cmax= max(abs(fc));
  87.   cmin= -cmin;
  88.   level= cmax/4.;
  89.   plf, fc, -y, -x, cmin=cmin, cmax=cmax;
  90.   animate, 1;
  91.   plf, fc, -y, -x, cmin=cmin, cmax=cmax;
  92.   plc, f, -y, -x, levs=0., marks=0, color="green", type="solid";
  93.   plc, f, -y, -x, levs=level, marks=0, color="black", type="dash";
  94.   plc, f, -y, -x, levs=-level, marks=0, color="green", type="dash";
  95.   fma;
  96.   pause, 1000;
  97.   tim_fill= array(0.0, 3);
  98.   timer, tim_fill;
  99.   tim_0= tim_fill;
  100.   for (i=1 ; i<=npass ; i++) {
  101.     lf= laplacian(f, y,x);
  102.     af= abs(f(2:-1,2:-1));
  103.     fdot+= lf*dt;
  104.     f(2:-1,2:-1)+= fdot*dt;
  105.     fc= f(zcen,zcen);
  106.     cmin= cmax= max(abs(fc));
  107.     cmin= -cmin;
  108.     plf, fc, -y, -x, cmin=cmin, cmax=cmax;
  109.     /* the 0 contour level is too noisy without some smoothing... */
  110.     fs= f(zcen,zcen)(pcen,pcen);
  111.     plc, fs, -y, -x, levs=0., marks=0, color="green", type="solid";
  112.     plc, f, -y, -x, levs=level, marks=0, color="black", type="dash";
  113.     plc, f, -y, -x, levs=-level, marks=0, color="green", type="dash";
  114.     fma;
  115.     timer, tim_fill;
  116.     if (tim_fill(3)-tim_0(3) > 60.) {
  117.       write, "aborting after one minute, "+pr1(i)+" frames";
  118.       break;
  119.     }
  120.   }
  121.   timer, tim_fill;
  122.   tim_fill -= tim_0;
  123.   pause, 1000;
  124.   animate, 0;
  125.   write,format="Rate for filled mesh is %f frames/(CPU sec), %f frames(wall sec)\n",
  126.     i/(tim_fill(1)+tim_fill(2)+1.0e-4), i/(tim_fill(3)+1.0e-4);
  127.  
  128.   /* roll the perspective movie */
  129.   persp: if (which && which!=2) goto shade;
  130.   f= f0;
  131.   limits, -1, 1, -1, 1;
  132.   pl3d,0, f, y, x;
  133.   animate, 1;
  134.   pl3d,0, f, y, x;
  135.   fma;
  136.   pause, 1000;
  137.   tim_3dw= array(0.0, 3);
  138.   timer, tim_3dw;
  139.   tim_0= tim_3dw;
  140.   for (i=1 ; i<=npass ; i++) {
  141.     lf= laplacian(f, y,x);
  142.     af= abs(f(2:-1,2:-1));
  143.     fdot+= lf*dt;
  144.     f(2:-1,2:-1)+= fdot*dt;
  145.     pl3d,0, f, y, x;
  146.     fma;
  147.     timer, tim_fill;
  148.     if (tim_fill(3)-tim_0(3) > 60.) {
  149.       write, "aborting after one minute, "+pr1(i)+" frames";
  150.       break;
  151.     }
  152.   }
  153.   timer, tim_3dw;
  154.   tim_3dw -= tim_0;
  155.   pause, 1000;
  156.   animate, 0;
  157.   write,format="Rate for 3D wire frame is %f frames/(CPU sec), %f frames(wall sec)\n",
  158.     i/(tim_3dw(1)+tim_3dw(2)+1.0e-4), i/(tim_3dw(3)+1.0e-4);
  159.  
  160.   /* roll the shaded movie */
  161.   shade: if (which && which!=3) return;
  162.   f= f0;
  163.   limits, -1, 1, -1, 1;
  164.   pl3d,1, f, y, x;
  165.   animate, 1;
  166.   pl3d,1, f, y, x;
  167.   fma;
  168.   pause, 1000;
  169.   tim_3dcol= array(0.0, 3);
  170.   timer, tim_3dcol;
  171.   tim_0= tim_3dcol;
  172.   for (i=1 ; i<=npass ; i++) {
  173.     lf= laplacian(f, y,x);
  174.     af= abs(f(2:-1,2:-1));
  175.     fdot+= lf*dt;
  176.     f(2:-1,2:-1)+= fdot*dt;
  177.     pl3d,1, f, y, x;
  178.     fma;
  179.     timer, tim_fill;
  180.     if (tim_fill(3)-tim_0(3) > 60.) {
  181.       write, "aborting after one minute, "+pr1(i)+" frames";
  182.       break;
  183.     }
  184.   }
  185.   timer, tim_3dcol;
  186.   tim_3dcol -= tim_0;
  187.   pause, 1000;
  188.   animate, 0;
  189.   pl3d,1, f, y, x;
  190.   fma;
  191.   limits;
  192.   write,format="Rate for 3D wire frame is %f frames/(CPU sec), %f frames(wall sec)\n",
  193.     i/(tim_3dcol(1)+tim_3dcol(2)+1.0e-4), i/(tim_3dcol(3)+1.0e-4);
  194. }
  195.  
  196. func laplacian(f, y,x)
  197. {
  198.   /* There are many ways to form the Laplacian as a finite difference.
  199.      This one is nice in Yorick.  */
  200.   /* Start with the two median vectors across each zone.  */
  201.   fdz= f(dif,zcen);
  202.   fzd= f(zcen,dif);
  203.   xdz= x(dif,zcen);
  204.   xzd= x(zcen,dif);
  205.   ydz= y(dif,zcen);
  206.   yzd= y(zcen,dif);
  207.  
  208.   /* Estimate the gradient at the center of each cell.  */
  209.   area= xdz*yzd - xzd*ydz;
  210.   gradfx= (fdz*yzd - fzd*ydz)/area;
  211.   gradfy= (xdz*fzd - xzd*fdz)/area;
  212.  
  213.   /* Now consider the mesh formed by the center points of the original.  */
  214.   x= x(zcen,zcen);
  215.   y= y(zcen,zcen);
  216.   xdz= x(dif,);
  217.   xzd= x(,dif);
  218.   ydz= y(dif,);
  219.   yzd= y(,dif);
  220.   area= xdz(,zcen)*yzd(zcen,) - xzd(zcen,)*ydz(,zcen);
  221.  
  222.   return ((xdz*gradfy(zcen,)-ydz*gradfx(zcen,))(,dif) +
  223.       (yzd*gradfx(,zcen)-xzd*gradfy(,zcen))(dif,)) / area;
  224. }
  225.  
  226. func pl3d(shading, z, y, x)
  227. {
  228.   /* rotate so that (zp,yp) are screen (y,x) */
  229.   /* These orientations are cunningly chosen so that the painter's
  230.      algorithm correctly draws hidden surfaces first -- see help, plf
  231.      for a description of the order cells are drawn by plf.  */
  232.   theta= 30. * pi/180.;  /* angle of viewer above drumhead */
  233.   phi= 120. * pi/180;
  234.  
  235.   ct= cos(phi);
  236.   st= sin(phi);
  237.   yp= y*ct - x*st;
  238.   xp= x*ct + y*st;
  239.  
  240.   ct= cos(theta);
  241.   st= sin(theta);
  242.   zp= z*ct - xp*st;
  243.   xp= xp*ct + z*st;
  244.  
  245.   if (!shading) {
  246.     color= [];
  247.     edges= 1;
  248.   } else {
  249.     /* compute the two median vectors for each cell */
  250.     m0x= xp(dif,zcen);
  251.     m0y= yp(dif,zcen);
  252.     m0z= zp(dif,zcen);
  253.     m1x= xp(zcen,dif);
  254.     m1y= yp(zcen,dif);
  255.     m1z= zp(zcen,dif);
  256.     /* define the normal vector to be their cross product */
  257.     nx= m0y*m1z - m0z*m1y;
  258.     ny= m0z*m1x - m0x*m1z;
  259.     nz= m0x*m1y - m0y*m1x;
  260.     n= abs(nx, ny, nz);
  261.     nx/= n;
  262.     ny/= n;
  263.     nz/= n;
  264.     color= bytscl(nx, cmin=0.0, cmax=1.0);
  265.     edges= 0;
  266.   }
  267.  
  268.   plf, color, zp, yp, edges=edges;
  269. }
  270.